knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.height = 3, fig.width = 12)
library(terra)
library(oharac)
library(tidyverse)
library(here)
source(here('common_fxns.R'))With impact maps generated for all stressors, means calculated both unweighted across species and FV-weighted across functional entities, we can add the impact layers to generate a comprehensive cumulative human impact map.
See scripts 2, 3a, and 3b for vulnerability mapping scripts, and scripts 4a-4d for impact calculation scripts.
Sum impacts across all stressor/vulnerability combos.
imp1_fs <- list.files(here_anx('_output/impact_maps/impact_maps_species'),
pattern = 'mean',
full.names = TRUE)
imp1_stack <- rast(imp1_fs) %>%
setNames(str_remove_all(basename(imp1_fs), '.+_spp_|_x_.+|_mean.tif'))ocean_a_rast <- rast(here('_spatial/ocean_area_mol.tif'))
### Calculate cumulative impacts
chi_spp_raw <- sum(imp1_stack, na.rm = TRUE)
chi_spp <- chi_spp_raw %>%
mask(ocean_a_rast) %>%
setNames('chi_spp')
writeRaster(chi_spp, here('_output/cumulative_impact_maps/chi_species.tif'),
overwrite = TRUE)### plot with ggplot
chi_spp_plot <- here('figs/chi_species.png')
chi_spp_df <- chi_spp %>%
as.data.frame(xy = TRUE) %>%
filter(!is.na(chi_spp))
fill_lbls <- c(0, .1, .25, .5, 1, round(minmax(chi_spp)[2], 2))
xfm <- function(x) {x^.25}
p <- ggplot(chi_spp_df, aes(x, y, fill = xfm(chi_spp))) +
geom_raster() +
scale_fill_viridis_c(labels = fill_lbls, breaks = xfm(fill_lbls)) +
coord_sf() +
theme_void() +
labs(title = 'CHI by species mean vulnerability',
fill = 'CHI')
ggsave(chi_spp_plot, height = 4, width = 8, dpi = 300)
knitr::include_graphics(chi_spp_plot)Sum impacts across all stressor/vulnerability combos.
imp2_fs <- list.files(here_anx('_output/impact_maps/impact_maps_funct_entity'),
pattern = 'mean',
full.names = TRUE)
imp2_stack <- rast(imp2_fs) %>%
setNames(str_remove_all(basename(imp2_fs), '.+_fe_|_x_.+|_mean.tif'))ocean_a_rast <- rast(here('_spatial/ocean_area_mol.tif'))
### Calculate cumulative impacts
chi_fe_raw <- sum(imp2_stack, na.rm = TRUE)
chi_fe <- chi_fe_raw %>%
mask(ocean_a_rast) %>%
setNames('chi_fe')
writeRaster(chi_fe, here('_output/cumulative_impact_maps/chi_funct_entity.tif'),
overwrite = TRUE)### plot with ggplot
chi_fe_plot <- here('figs/chi_funct_entity.png')
chi_fe_df <- chi_fe %>%
as.data.frame(xy = TRUE) %>%
filter(!is.na(chi_fe))
fill_lbls <- c(0, .1, .25, .5, 1, round(minmax(chi_fe)[2], 2))
p <- ggplot(chi_fe_df, aes(x, y, fill = xfm(chi_fe))) +
geom_raster() +
scale_fill_viridis_c(labels = fill_lbls, breaks = xfm(fill_lbls)) +
coord_sf() +
theme_void() +
labs(title = 'CHI by functional entity vulnerability',
fill = 'CHI')
ggsave(chi_fe_plot, height = 4, width = 8, dpi = 300)
knitr::include_graphics(chi_fe_plot)Difference in CHI calculated from FE vulnerability relative to species vulnerability (using only spp included in functional entities): \[\frac{CHI_{FE} - CHI_{spp}}{CHI_{spp}} \times 100\%\]
r <- ((chi_fe - chi_spp) / chi_spp * 100) %>%
setNames('layer')
r_vals <- values(r); r_vals <- r_vals[!is.na(r_vals)]
r_zlim <- max(abs(quantile(r_vals, .001)), abs(quantile(r_vals, .999)))
values(r)[values(r) > r_zlim] <- r_zlim
values(r)[values(r) < -r_zlim] <- -r_zlim
### set up thresholded rasts as well
r1 <- r2 <- r
r1[chi_fe < .10 & chi_spp < .10] <- NA
r1[r1 > r_zlim] <- r_zlim
r1[r1 < -r_zlim] <- -r_zlim
r2[chi_fe < .25 & chi_spp < .25] <- NA
r2[r2 > r_zlim] <- r_zlim
r2[r2 < -r_zlim] <- -r_zlim
### Difference plots
chi_diff_plot <- here('figs/chi_fe_diff_spp.png')
chi_diff1_plot <- here('figs/chi_fe_diff_spp_0.10.png')
chi_diff2_plot <- here('figs/chi_fe_diff_spp_0.25.png')
chi_diff_df <- as.data.frame(r, xy = TRUE) %>%
filter(!is.na(layer))
map_cols <- hcl.colors(palette = 'Red-Green', n = 15, rev = TRUE)
p <- ggplot(chi_diff_df, aes(x, y, fill = layer)) +
geom_raster() +
scale_fill_gradientn(colors = map_cols,
limits = c(-r_zlim, r_zlim)) +
coord_sf() +
theme_void() +
labs(title = 'CHI difference between FE and species',
fill = '%∆CHI')
ggsave(chi_diff_plot, height = 4, width = 8, dpi = 300)
### Difference cropped by .10 threshold
chi_diff1_df <- as.data.frame(r1, xy = TRUE) %>%
filter(!is.na(layer))
p <- ggplot(chi_diff1_df, aes(x, y, fill = layer)) +
geom_raster() +
scale_fill_gradientn(colors = map_cols,
limits = c(-r_zlim, r_zlim)) +
coord_sf() +
theme_void() +
labs(title = 'CHI difference between FE and species (0.10 threshold)',
fill = '%∆CHI')
ggsave(chi_diff1_plot, height = 4, width = 8, dpi = 300)
### Difference cropped by .25 threshold
chi_diff2_df <- as.data.frame(r2, xy = TRUE) %>%
filter(!is.na(layer))
p <- ggplot(chi_diff2_df, aes(x, y, fill = layer)) +
geom_raster() +
scale_fill_gradientn(colors = map_cols,
limits = c(-r_zlim, r_zlim)) +
coord_sf() +
theme_void() +
labs(title = 'CHI difference between FE and species (0.25 threshold)',
fill = '%∆CHI')
ggsave(chi_diff2_plot, height = 4, width = 8, dpi = 300)
knitr::include_graphics(chi_diff_plot)In and outside of neritic zone.
neritic_rast <- rast(here('_spatial/bathy_mol_neritic.tif'))
mean_str_spp <- sapply(imp1_stack,
FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
setNames(names(imp1_stack))
mean_str_spp_coastal <- sapply(imp1_stack %>% mask(neritic_rast),
FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
setNames(names(imp1_stack))
mean_str_fe <- sapply(imp2_stack,
FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
setNames(names(imp2_stack))
mean_str_fe_coastal <- sapply(imp2_stack %>% mask(neritic_rast),
FUN = function(r) mean(values(r), na.rm = TRUE)) %>%
setNames(names(imp2_stack))
mean_str_df <- data.frame(stressor = names(imp1_stack)) %>%
mutate(mean_impact_fe = mean_str_fe[stressor],
mean_impact_fe_nonzero = mean_str_fe_coastal[stressor],
mean_impact_spp = mean_str_spp[stressor],
mean_impact_spp_nonzero = mean_str_spp_coastal[stressor]) %>%
arrange(desc(mean_impact_fe))
knitr::kable(mean_str_df)| stressor | mean_impact_fe | mean_impact_fe_nonzero | mean_impact_spp | mean_impact_spp_nonzero |
|---|---|---|---|---|
| sst_extremes | 0.0338281 | 0.0488863 | 0.0351647 | 0.0491788 |
| ocean_acidification | 0.0330115 | 0.0567399 | 0.0244016 | 0.0547503 |
| sst_rise | 0.0270853 | 0.0242258 | 0.0228888 | 0.0234100 |
| uv_radiation | 0.0152301 | 0.0306834 | 0.0172201 | 0.0303766 |
| direct_human | 0.0092048 | 0.0092048 | 0.0087021 | 0.0087021 |
| biomass_removal | 0.0078020 | 0.0241021 | 0.0096138 | 0.0273267 |
| sea_level_rise | 0.0057294 | 0.0057294 | 0.0058910 | 0.0058910 |
| bycatch | 0.0040566 | 0.0155488 | 0.0042723 | 0.0151240 |
| shipping_large | 0.0023091 | 0.0108484 | 0.0025661 | 0.0121029 |
| light | 0.0010080 | 0.0100309 | 0.0010055 | 0.0099937 |
| nutrient | 0.0008922 | 0.0090952 | 0.0008932 | 0.0091172 |
| fishing_benthic_dest | 0.0001761 | 0.0016030 | 0.0001688 | 0.0015191 |
| benth_str | 0.0000046 | 0.0000449 | 0.0000043 | 0.0000415 |